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We discuss the effects of quantum fluctuations on the properties of vortex lattices in rapidly rotating ultracold 
atomic gases. We develop a variational method that goes beyond the Bogoliubov theory by including the effects 
of interactions between the quasiparticle excitations. These interactions are found to have significant quantitative 
effects on physical properties even at relatively large filling factors. We use our theory to predict the expected 
experimental signatures of quantum fluctuations of vortices, and to assess the competition of the triangular 
vortex lattice phase with other phases in finite-sized systems. 
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I. INTRODUCTION 

The achievement of Bose-Einstein condensation in dilute 
atomic gases has allowed studies of the interesting collec- 
tive behaviour of these quantum fluids with unprecedented 
detail and in previously inaccessible parameter regimes J2]. 
The interatomic interactions, although weak, cause the atomic 
Bose-Einstein condensates to exhibit superfluidity and the as- 
sociated phenomena familiar from studies of superfluid he- 
lium. Among the most striking of these features is the ap- 
pearance of quantized vortices when the superfluid is forced 
to rotate |B|3|]. Experiments on rapidly rotating atomic gases 
have shown beautiful images of vortex lattices flj-la], and have 
allowed detailed studies of their collective properties 10, [H]. 

Theoretical studies have pointed to the possibility of using 
ultracold atomic gases to explore novel regimes of vortex den- 
sity where quantum fluctuations of the vortices become im- 
portant The parameter controlling the degree of quantum 
fluctuations is the filling factor, v = N/N v , where N and N v 
denote the number of bosons and vortices respectively ifioll . 
Quantum fluctuations are small for v — s- oo, but increase with 
decreasing v. In particular, for a large uniform system, it has 
been shown that these fluctuations drive a phase transition 
out of the triangular vortex lattice to a set of strongly corre- 
lated bosonic quantum liquids (analogous to fractional quan- 
tum Hall states) when the filling factor is reduced below a crit- 
ical value v c IU0I1 . Estimates of v c vary between u c ~ 2 — 6 
from exact diagonalization results | llTj|4l2ll and v c ~ 8 — 14 
from a Lindemann criterion ITol [l3 Jl4ll . 

While the prediction of the precise transition point is a very 
difficult theoretical task, the analytical theories of Refs. 11131 — 
[Till do also provide important information concerning the ef- 
fects of quantum fluctuations within the vortex lattice phase. 
These results are based on the application of mean-field the- 
ory, with quantum fluctuations described by the Bogoliubov 
approximation. This is expected to be accurate in the limit of 
large filling factors, v 3> 1 ifioll . However, for smaller values 
of v, one expects corrections to the results of Bogoliubov the- 
ory, arising from the effects of quartic terms in the fluctuation 
expansion, or equivalently from effective interactions between 
the Bogoliubov quasiparticles. In this paper, we analyse the 
effects of these corrections to the Bogoliubov theory. We in- 
troduce a variational wavefunction for the vortex lattice of a 
rapidly rotating atomic Bose gas, which recovers the Bogoli- 



ubov theory at v — > oo but which incorporates the effects of 
interactions between the Bogoliubov quasiparticles. We show 
that the variational theory accurately reproduces the exact di- 
agonalization results for v > 5 in the small systems for which 
these results are available. We then apply the theory to large 
systems, representative of the thermodynamic limit. We find 
that the quasiparticle interactions lead to large quantitative 
corrections to the predictions of Bogoliubov theory even for 
large filling factors deep in the vortex lattice regime. We dis- 
cuss how effects of quantum fluctuations may be observed in 
experiments measuring the condensate fraction, particle den- 
sity at the vortex core or the shear modulus of the vortex lat- 
tice. Finally, we discuss the implications of our results for the 
competition with other phases on finite-sized systems. 

The paper is organised as follows: In Sec. [U]we review the 
basic model for our system and derive the form of the mi- 
croscopic Hamiltonian in the rapid-rotation limit. We then 
apply Bogoliubov theory in Sec. [Ill] obtaining the spectrum 
of low-energy excitations for the triangular lattice, and using 
this to derive the effects of quantum fluctuations on various 
physical observables. In Sec. IIV1 we present a variational 
approach which allows studies beyond the Bogoliubov limit 
of very small condensate depletion. In Sec |VJ we compare 
the results of this variational theory with small system exact 
diagonalization results, and describe its predictions both for 
(thermodynamically) large systems and for finite-sized sys- 
tems. Sec|VT]contains a summary of our results. 



H. THEORETICAL MODEL 

A system of N interacting bosons of mass M in a parabolic 
trap that is rotating with frequency SI about the z-axis can be 
described by the following rotating-frame Hamiltonian @]: 
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where Vt r and 2 are the radial and axial trap frequencies 
respectively and V(r) = g3uS( r ) i s a contact interaction 
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of strength g^. The Hamiltonian is equivalent to a sys- 
tem of charge-g bosons in an effective magnetic field B = 
(2MSl/q)z. We will be working under the assumptions of 
weak interactions and quasi 2-dimensionality where the in- 
equalities ng -C Ml and ng -c Ml z are satisfied (n is the 

2D boson density and g 



93D 



2-Kh 



is the 2D contact in- 



teraction parameter) 11711 . We set fl r = O so that any resid- 
ual radial confinement can be neglected. Treating the interac- 
tion term as a perturbation, it follows that the energy eigen- 
states of H are Landau levels with eigenvalues given by 
E n = 2Ml(n + 1/2), and that all the particles are confined 
to the lowest Landau level (LLL). We shall refer to the above 
theoretical set-up as the 2D LLL regime. The LLL basis states 
take the following form: 
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where I = yj K/2MQ, is the effective magnetic length and 
j = 0, 1, 2, . . . indexes the states. 

Within mean-field theory, the condensate wavefunction, 
'i'(r), is a superposition of these basis states and is determined 
solely by the positions of its zeros IU8I1 . These correspond to 
coordinates of vortex centres. Under the requirement that the 
entire system mimics rigid body motion the number of vor- 
tices in the system must be given by 



N v = A/(2vr/ 2 ), 



(3) 



where A is the condensate area. 

Studying the structure of the mean-field ground state will 
help us to choose an appropriate basis in which to treat quan- 
tum fluctuations. In the 2D LLL regime we are left to consider 
only the interaction term in the Hamiltonian. In the mean-field 
approximation the interaction term reduces to 
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|*(r)| 4 d 2 r, 
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and assumes a minimum when the zeros of ^(r) lie on a tri- 
angular lattice |9|]. Thus, the mean-field ground state is a tri- 
angular lattice of vortices with a unit cell area of 2tt1 2 (a con- 
sequence of Eq. [3]l. 

As our basis, we will use LLL magnetic Bloch states which 
are related to the mean-field ground state by a translation. We 
review their construction, outlined in Ref. lfl9ll . One starts 
with the most localized symmetric LLL wavefunction 



co(r) = $j=o(r) 
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and translates it to the sites of a 2D Bravais lattice, to obtain 
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where a! and a 2 are the Bravais lattice vectors, r m = 7774 a! + 
m2a 2 is the translation vector from the origin to a given lat- 
tice site, and T rm = cxp [— ^r m ■ (p + milz, x r)] is the 



corresponding translation operator rf20ll . From Eq. [3] it fol- 
lows that one quantum flux passes through the unit cell area 
|ai x &2\, and therefore, the magnetic translation operators 
commute with the Hamiltonian and each other, generating the 
elements of the Magnetic Translation Group |20Q . Following 
Ref. lfl9ll . we can now construct orthonormal Bloch functions 
out of linear combinations of c m (r) 



*k(r) 
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where £(k) normalises the Bloch function to unity over the 
system area, and is given by 



C(k) = ^(-l)™i™2 e -r m 2 /4Z 2 e - 4 k.r^ 
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The Bloch states are simultaneous eigenstates of the Hamilto- 
nian and the translation operators. In a finite system of area 
A, there are N v (Eq. [5]) such states, labelled by momenta k, 
which belong to the first Brillouin zone (BZ) of the 2D Bravais 
lattice and satisfy periodic boundary conditions with respect 
to magnetic translations. We shall study a system of vortices 
defined on a rectangular region with dimensions L x , L y , and 
with periodic boundary conditions. 

In the magnetic Bloch basis the contact interaction takes the 
following form: 

Hl= M ( k 'q-q''q") (S q+q',q"+Gfck +q fok +q ^k+q''&k, 
k,q,q',q",G 

(9) 

where all momenta belong to the first BZ and G are the recip- 
rocal lattice vectors. The Kronecker delta imposes conserva- 
tion of momentum modulo a reciprocal lattice vector and the 
interaction matrix element, obtained in Ref. (2j]|, is given by 

M(k,q,q',q") = 9 - J d 2 rvI/* +q (r)vI/* +q ,(r)vI/ k+q „(r)vI- k (r) 
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The mean-field ground state corresponds to all N particles 
condensed into the k = Bloch state. It has energy E = 

N 2 M(0, 0, 0, 0) =/3 A ^-,where^ A ~ 1.1596. 

One can readily extend this approach to calculate the shear 
modulus within mean-field theory. The energy change per unit 
area of a vortex lattice, sheared by an infinitesimal angle 8, is 
given by C-20 2 , where C2 is its shear modulus II 1 511 . The shear 
modulus can be calculated in the mean-field approximation as 
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where the 6 dependence of M(0, 0, 0, 0) enters through the 
sheared lattice vectors: 



ai = oa (1,0) 



a 2 = ^(1 + 0^3, V3). 



(12) 



a a = (^)il is the lattice constant of the triangular vortex 
lattice. The shear modulus is relevant for the distortion of 
the lattice by thermal fluctuations, and ultimately the thermal 
melting transition. (Though non-linearities from thermal fluc- 
tuations will become important.) It also determines the fre- 
quencies of the long-wavelength sound modes in the hydro- 
dynamic regime. This limit applies when the scattering rate 
Tk of Bogoliubov modes is large compared to their energy 
6k- Matveenko et al. lfl6ll have shown that at zero temper- 
ature, ftr(e(k), 0)/e(k) ~ 0.065/^ so the system is always 
lightly damped at filling factors where the triangular lattice is 
predicted to exist (i.e. v > 7). However, at non-zero temper- 
ature, excitations with e(k) < e c = are overdamped (see 



11611 ). Thus, for small non-zero temperatures, one can expect 
these modes to be in the hydrodynamic regime, and have fre- 
quencies set by the shear modulus. 

It is useful to note that the zeros of the Bloch function 
* k (r) are located at JlH 



1 , 
+ 5 (ai 



a 2 ) + l 2 z x k. 
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From this, we see that the zeros map out a triangular lattice 
which is simply a translation of the mean-field k = lattice. 
Adding quantum fluctuations to the mean-field theory leads to 
excitations of neighbouring k-states, which correspond to vor- 
tex lattices that are offset with respect to the macroscopically 
occupied one. We can say that quantum fluctuations "smear" 
the original mean-field lattice. 



III. BOGOLIUBOV APPROXIMATION 

At high filling factors, the condensate is weakly depleted, 
and we can approximate the Hamiltonian by a Bogoliubov ex- 
pansion, which keeps terms up to second order in (N — No), 
where No is the macroscopic occupation number of the con- 
densate. The condensate is a coherent state, i.e. an eigenstate 
of bo. Its creation and annihilation operators can be treated as 
c numbers given by 



b\ = bo = v^Vo= N-Y,(b& 
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The Bogoliubov expansion of Hi takes the following explicit 
form: 

Hi = AI(0,0,0,0)N 2 

+ o 1^ \ b v b -^) \ |A q | e -^ f„ I [ fit _ 
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where 

|A q |e 1 ^ = 2M(0,q,-q,0), 

|A q |e-** = 2M(q,-q,-q,-2q), 

£ q = 2 [M (0, 0, q, q) + M (0, q, 0, q) - M(0, 0, 0, 0)] . 

(16) 

In this approximation the Hamiltonian is quadratic and can be 
diagonalized via a Bogoliubov transformation of the form: 
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Hi = N 2 M(0, 0, 0, 0) 



q#0 



a-l A ql 2 K a q + ^)-^q 



(18) 

where tanh26L = The number of particles outside the 
condensate is given by 



N-No = 5Z(6 q 6 q ) = 5>inh 2 

q^O q#0 
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(19) 



which is assumed small compared to N for self-consistency of 
the mean field approach. However, in the vicinity of q = 0, 
one finds that sinh 2 q ~ l/(ql) 2 . Using £ q ->• ^ / d 2 q, 
one finds that the depleted fraction diverges with system 
size as (N - N )/N ~ ln(/V v )/iA 

The excitation spectrum of the triangular lattice, given by 
Eq.[l8] agrees with previous calculations lfl3i[l6ll . We extend 
these works on the Bogoliubov approximation by evaluating 
also the zero-point energy, the particle density at the vortex 
core, and the shear modulus. To do so, we calculate the matrix 
elements, £ q and |A q |, by evaluating the sums in Eq. [TOlfor 
the triangular lattice. (See the appendix for an outline of this 
calculation.) 

We find that including the zero-point fluctuations within 
Bogoliubov theory leads to the ground state energy 



(20) 



where 7a — —0.5036. Zero-point fluctuations give rise to 
a correction that increases as we move away from the mean- 
field regime and is inversely proportional to the filling factor. 

Condensate depletion, which takes place beyond mean- 
field, gives a non-zero expectation value of the particle density 
at the centre of the vortex: 
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where r r 



= |(ai + a 2 ) gives the position of one of the ze- 
ros of the mean-field k = wavefunction. In the Bogoliubov 
approximation, we find 
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f c ~ 5. As we shall discuss below, the inclusion of effects 
beyond the Bogoliubov theory have an even more dramatic 
effect. 

The zero-point fluctuations also give a significant correc- 
tion to the shear modulus, arising from the variation of the 
zero-point energy with the shear angle. Under the Bogoliubov 
approximation we find 
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where the 8 dependence enters through the sheared lattice vec- 
tors (see Eq.[T2l. 



FIG. 1: Vortex core particle density as a function of the inverse 
filling factor, for the full Bogoliubov theory (solid line), the long- 
wavelength theory of Ref. Ill3ll (dashed line), and for our variational 
wavefunction in Eq. [28] for N v = 90 (diamonds). The relation 
given in the appendix was used to map from the mean-squared vortex 
core displacement (Sr 2 ) calculated in Ref. lfl3ll to vortex core den- 
sity. n c (dotted line) is the critical core density at which the lattice 
should "melt" according to a Lindemann criterion with parameter 
(Sr 2 ) = 0.145/ 2 . 



The vortex core density can be used to estimate v c , the crit- 
ical filling factor at which the vortex lattice melts, by relating 
it to spatial fluctuations of the mean-field lattice. As was ex- 
plained in the preceding section, condensate depletion leads 
to excitation of k ^ 0-states, and this can be viewed as the 
fluctuating mean-field lattice. The fluctuations lead to a non- 
vanishing particle density at the core of each vortex. Assum- 
ing that the fluctuations are described by a Gaussian distribu- 
tion we can write down the average vortex core density as a 
function of the variance of this distribution, (Sr 2 ), where Sr 
is the displacement of the vortex core from its mean-field po- 
sition (see the appendix for this relation). For \8r\ <C I the 
density is a function of the variance only and does not de- 
pend on the particular distribution of the fluctuations. One 
can estimate the melting point of the lattice using a suit- 
able Lindemann criterion |22y i.e. the lattice is said to melt 
when (8r 2 )/l 2 reaches a particular threshold value. We fol- 
low Refs. IS [3] and take this to be ~ 0.145. We also cal- 
culate the vortex core density from the estimate of Ref. lfl3ll 
for (Sr 2 ), obtained through a path-integral formalism. Fig.Q] 
compares the core density in this long-wavelength limit with 
the result that we find within the Bogoliubov approximation. 
The Bogoliubov theory is more accurate because it uses the 
energy dispersion across the entire Brillouin zone, whereas the 
long-wavelength approximation of Ref. [130 uses a quadratic 
approximation for the energy dispersion which neglects the 
contribution of higher order k terms to (Sr 2 ). The inclusion 
of the full dependence across the Brillouin zone leads to a 
change in the melting transition predicted by the Lindemann 
criterion with parameter (5r 2 )/l 2 = 0.145 from v c ~ 9 to 



IV. BEYOND THE BOGOLIUBOV APPROXIMATION 

The Bogoliubov expansion of the interaction Hamiltonian 
(Eq. |9]l neglects those scattering processes which do not in- 
volve condensed particles, i.e. terms of the form & k & q &q'fr q ", 
and is only valid for small condensate depletions when v 1. 
Previous studies, summarized in Section [Till have overlooked 
the contribution of this quartic term. As will be shown, the 
contribution of the quartic term to the energy of the system 
can be very significant even at large filling factors. In fact, the 
larger the system the higher the filling factor at which we ex- 
pect the quartic term to be a significant contribution to the Bo- 
goliubov approximation. This is a result of the diverging de- 
pletion in the vicinity of k = 0, where sinh 2 (9k ~ N v within 
Bogoliubov theory. A scaling argument shows that this di- 
vergence causes the contribution of the quartic terms to the 
total energy to be of order N v /v times the contribution of 
the quadratic terms. Hence, the Bogoliubov approximation 
is quantitatively valid, and the quartic term can be neglected, 
only when N v /v <C 1. 

The most straightforward way to include the quartic term 
in energy calculations is through a variational analysis. We 
take a trial wavefunction of the form obtained within the Bo- 
goliubov approximation (retaining up to quadratic terms), pa- 
rameterised by a set of variational variables. We then find 
the values of these parameters that optimise the full energy, 
including the quartic interactions. The Bogoliubov wavefunc- 
tion can be expressed in terms of a depletion operator acting 
on a simple condensate. The operator excites and de-excites 
pairs of particles with a vanishing total momentum. The exci- 
tations are weighted by a set of real variational parameters, #k 
and 0k, indexed by the BZ momenta. 



Atrial) = U\0), 
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where |0) is the condensate of all particles in the k = state, 
and U a unitary depletion operator given by 



U = cxp 
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By making use of the following identities: we can express (Atrial I -Hi I Atrial) as an explicit function of the 

parameters 8^ and tp^: 

U^bqU = &qcosri0 q - bl^o smh6» q , 
U^b\U = b\ cosh 6» q - b.qe-^ sinh 6> q , (26) 



(*trial|^i| Atrial) = J2 {2 sinh 2 k sinh 2 dyif (k, 0, -k + q, -k + q) 

k^0,q^0 

+M(0, 0,0, 0)TV 2 + sinh 6» q cosh 6> q sinh 6> k cosh 6» k e i(0k ~ 0q) M(k, k + q, q k, 2k) | 

+7V ^] [4 sinh 2 6» q A/(0,q,0,q) -sinh 6» q cosh 6» q e^M (q, -q, -q, -2q) + e^"A/(0, q, -q, 0) j, (27) 
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where TVo can be eliminated using Eq. [19] 

The Bogoliubov trial wavefunction does not have a definite 
particle number. The finite variance of TV leads to a significant 
overestimate of the energy at low filling factors and in small 
systems, i.e. whenever TV is not much greater than 1. We 
can see this from how the energy scales with TV in the mean- 
field regime: ^j-N 2 . It follows that E/N 2 is overestimated 



by 
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To address this problem we could project 



the Bogoliubov trial wavefunction onto a state with a definite 
particle number and compute the new expectation value of H\ . 
We instead propose to use, right from the beginning, a trial 
wavefunction with a definite particle number of the following 
form, in line with work by Girardeau and Arnowitt l23ll : 
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where fii = bo(blbo)~ 2 and |TV) is a condensate of TV parti- 
cles in the k = state. The depletion operator is now particle- 
conserving. The expectation value of H\ for this trial wave- 
function is expressed in terms of the variational parameters 
in the appendix. The optimal values of the parameters were 
found numerically for each system that we studied. 



V. RESULTS AND COMPARISONS WITH OTHER STATES 

A. Triangular Vortex Lattice 

In order to test our variational method we have calculated 
the energies of small systems and have compared them with 
the data available from exact diagonalisation (ED) studies of 
Ref. [10]. Variational energies become exact when v — > 00 so 
should match the ED results in that limit. Fig. [2] presents 
our variational bounds, obtained from both the indefinite- 
number and the definite-number trial wavefunctions, and com- 
pares them against the ED energies. Clearly, the state with 
definite particle number provides a much better variational 
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FIG. 2: Comparison of our results with ED data for TV V = 6 and 
L y /L x — l/\/3. As expected and explained in section ITVl the in- 
definite particle number trial wavefunction overestimates the energy 
and this overestimate is larger at low filling factors. In all the other 
figures our variational results were obtained for the definite number 
trial wavefunction in Eg. 1281 



bound. The energies of the definite-number state are accu- 
rate to within 1%, down to v ~ 5. (For filling factors v < 5 
we expect the groundstate to involve significant quantum fluc- 
tuations, being close to or in the regime of strong correlated 
quantum Hall states, and the variational wavefunction for the 
vortex lattice to become inaccurate. Still, we present our re- 
sults for the range < 1/v < 1 for completeness.) 

We can now take advantage of this success and study the 
behaviour of the condensate for large systems in the thermo- 
dynamic limit, well outside the range of the available ED data. 

Fig. [3] shows the results of variational calculations for the 
largest system we have studied, of TV V = 90 vortices. The 
lowest energy variational state is the triangular lattice (shown 
by triangular symbols in Fig. [3]). At large filling factors, 
1/v 0, the variational results match with the Bogoliubov 
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FIG. 3: The results for the largest studied system, i.e. for N v = 90 
(Ly/L x = |v3 for the triangular lattice and L y /L x — 9/10 for 
the square lattice). The variational energies were calculated using 
the definite-number trial wavefunction. 



FIG. 4: Condensate depletion in a range of system sizes. The sym- 
bols correspond to our variational results for the definite-number 
state and the solid lines show the depletion calculated in the Bogoli- 
ubov approximation. 



theory. However, there are significant deviations from the Bo- 
goliubov result at low filling factors. For v ~ 10 the deviation 
is as large as the difference in energy we find between trian- 
gular and square lattices, showing that these non-linearities 
are quantitatively important in determining the nature of the 
low-energy phase. 

We have also computed the condensate depletion, in a range 
of system sizes, as a function of the filling factor v (see Fig.@]i. 
The inclusion of the quartic term quenches condensate deple- 
tion significantly. Already, at v ~ 10, we see a 15% difference 
in condensate depletion for iV v = 90 when it is included. The 
condensate depletion, although substantially reduced from the 
prediction of Bogoliubov theory, is still significant. An ex- 
perimental observation of this reduction from 100% conden- 
sate fraction - even in the limit of vanishing temperature - 
would be a clear indication of the effects of quantum fluctua- 
tions in. 

As discussed in Sec|Illl condensate depletion leads to a non- 
vanishing vortex core density. Fig. [T] shows the vortex core 
density for iV v = 90, obtained from the variational estimates 
of the occupation numbers of k ^ states. The contribution 
of the quartic term is very significant and cannot be neglected 
at filling factors of the order of iV v ~ 90, where Bogoliubov 
is no longer a valid approximation (see the scaling argument 
in Section. HV]). 

It is particularly striking that, for filling factors as small as 
v = 1, the core density does not reach the critical value n c 
obtained from the Lindemann melting criterion with parame- 
ter (<5r 2 ) = 0.145Z 2 . The application of this numerical value 
for the Lindemann criterion is highly questionable in this case 
of quantum melting of a vortex lattice. Since it is known from 
numerical exact diagonalization studies that the lattice melts 
with v c ~ 2 — 6 IHM12I1 . our results for the core density sug- 
gest that the Lindemann parameter should be significantly re- 
duced. We do not believe it is helpful to quote a number here, 
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FIG. 5: Shear modulus (C2) of the triangular vortex lattice as a func- 
tion of the filling factor for iV v = 90 



since the variational approach (based on expanding around the 
triangular vortex lattice) may be inaccurate in the vicinity of 
the melting transition. Nevertheless, note that there is already 
a very clear departure from the Bogoliubov results for large 
v ~ N v . Our results show that the effects of quantum fluctua- 
tions in the positions of vortices are very much suppressed as 
compared to the expectations of Bogoliubov theory. 

Including quantum fluctuations beyond the Bogoliubov ap- 
proximation also leads to very significant corrections to the 
shear modulus of the vortex lattice. Fig. [5] shows our numer- 
ical values of the shear modulus for iV v = 90 and a range of 
filling factors, and compares them against the Bogoliubov ap- 
proximation. The shear modulus is not suppressed as much 
as predicted by Bogoliubov theory in Eq. [23] The shear mod- 
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FIG. 6: Energy of the smectic state (circles) compared against the en- 
ergy of the triangular vortex lattice (triangles) for N v = 8, L y /L x — 
V§/4 and for N v = 24, L y /L x = 4\/3/9 



ulus can be experimentally measured from the frequencies of 
the Tkachenko sound modes. Following the continuum the- 
ory lfl5l l25ll of the Tkachenko modes in a rapidly rotating 
BEC, the eigenfrequencies are given by 

/ 47rC9 

where 71 = 7.17 and 72 = 16.9 for the two lowest eigen- 
modes. The Tkachenko frequencies are proportional to \/C2- 
As described above, depending upon whether the oscillations 
are much faster or slower than the damping rates of Bogoli- 
ubov excitations, the system is in the collisionless or hydro- 
dynamic regime. In the hydrodynamic regime zero-point fluc- 
tuations need to be taken into account. Our results show a 
4% shift from the mean-field value of C2, at filling factors 
of ~ 10, giving an approximately 2% shift in the Tkachenko 
eigenfrequencies. Observing this small softening of the shear 
modulus below its mean-field value would be an indication of 
the role of quantum fluctuations in the vortex lattice phase. 



B. Competing Phases in Small Systems 

Finally, we turn to discuss the energetic stability of the tri- 
angular lattice phase compared with other competing phases 
in finite-sized systems. One possible phase the sytem might 
enter, as the filling factor or system size are reduced, is the 
smectic phase. In Ref. iflltl the smectic state was shown to 
be favoured over the triangular lattice in small systems. The 
smectic state corrresponds to N/N s independent condensates 
put into Landau-gauge states, which have the form of stripes 
aligned along the width of the system (7V S is the number of 
particles put into each state). We have compared our results 
for the vortex lattice phase with the energies of the stripe phase 



in finite-sized systems. Following Ref. ifTHl . the stripe separa- 
tion that optimises the energy is given by 0.9083q,a- For each 
system we studied, we chose the separation that is closest to 
this value and commensurate with the length of the system, 
L x . In the presence of a contact interaction, we only need to 
consider nearest-neighbour interactions between the stripes. 
This means that (Hi) /N 2 depends solely on the number of 
particles per stripe, N s , and the separation of the stripes, Ax. 
For N v = 4, 6, 8, where Ax = and there are two vortices 
per stripe (N v , s = 2), the energy of the smectic state becomes 

(H)/N 2 = 2^(1-1888 - 1.0746/Ay, (30) 
for N v = 24, where Ax = ^a A and iV v>s = 4, 

(Hi)/N 2 = ^-(1.1757 -0.9307/JVb), (31) 
and for iV v = 90, where Ax = j^oa and N V)S = 9, 

(Hi)/N 2 = ^-(1.1720 -0.9671/iV s ). (32) 

Our results agree qualitatively with those of Ref. fPUl . In 
the thermodynamic limit the triangular vortex lattice has lower 
energy than the smectic state down to small filling factors at 
which we expect to see fractional quantum Hall states iflfjll 
(FQHSs). We also find that the smectic phase is favoured for 
small systems and low filling factors, and that there is evi- 
dence for phase transitions between the smectic state and the 
triangular vortex lattice. However, there are important quanti- 
tative differences from the results of Ref. ifTHl . That work did 
not include quantum fluctuations, obtaining ^j-{N 2 — N) 
for the energy of the triangular vortex lattice. By including 
the Bogoliubov and quartic terms in our energy calculations, 
we revise the estimates of the transition points between the 
triangular and smectic states. For N v =4,6,8 (N VyS = 2) our 
results show a phase transition at v = 3 — 4. As expected, 
quantum fluctuations stabilise the triangular vortex lattice and 
our estimate is lower than v = 8 — 13, obtained following 
Ref. ifTlll . Indeed, we find that once the system size is as large 
as N v = 24 (see Fig. [5), there is no transition to a smectic 
state for any filling factors v > 1. 

Our studies show that another competing phase is the 
square vortex lattice. Like the smectic state, it can be favoured 
over the triangular lattice in small systems and at low filling 
factors. Following the same variational method as for the tri- 
angular lattice, we computed the matrix elements in Eq. [10] 
with the lattice vectors now describing a square lattice, and 
calculated variational energies using the definite-number state 
(Eq.|28li. We found an energy cross-over between the two vor- 
tex lattice phases for N v = 6 at v ~ 4. In the largest, approx- 
imately isotropic, system that we studied, i.e. N v — 90 (see 
Fig. |3), the triangular vortex lattice is lower in energy down 
to v ~ 1. Thus, our results show that, while quantum fluctua- 
tions reduce the energy gap between the square and triangular 
vortex lattice phases (and, as described above, also reduce the 
shear modulus of the triangular lattice), the triangular vortex 
lattice remains the lower energy state in the thermodynamic 
limit. 
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VI. SUMMARY 

We have developed a variational method to describe the 
quantum fluctuations of 2D vortex lattice phases of ultracold 
atomic gases in the lowest Landau level regime. We find very 
good agreement (within 1%) between the results of our vari- 
ational method and those of exact diagonalization studies at 
filling factors of 5 and above. Our theory includes the effects 
of interactions between the quantum fluctuations beyond the 
Bogoliubov description. These are found to have significant 
quantitative effects on physical properties. Our theory pre- 
dicts dramatically reduced condensate depletion as compared 
to Bogoliubov theory, and dramatic reduction in the extent to 
which the positions of the vortices fluctuate, as evidenced by 
the quantum contribution to the particle density at the vortex 
cores. Measurements of the condensate fraction, of the shear 
modulus, or of the particle density at the vortex cores could 
provide signatures of the quantum fluctuations of vortices. 
Our results provide clear predictions for the sizes of these 
effects, which differ markedly from the predictions of previ- 
ous theories even at relatively large filling factors. Finally, 
we used our theory to assess the competition of the triangu- 
lar vortex lattice phase with the smectic and square phases in 
finite-size systems. 



ansatz that vortex core positions are described by the follow- 
ing Gaussian probability distribution: 



p(6r) 



1 



7r((5r 2 



-5r 2 /(5r 2 ) 



(Bl) 



where Sr is the displacement of the vortex core from its mean- 
field position and (<5r 2 ) = ^l 2 is the variance of that displace- 
ment. It turns out, however, that the relation we are deriving 
is independent of the above ansatz when the vortex core fluc- 
tuations are small \Sr\ -C I, which is the case right up to the 
melting point given by the Lindemann criterion that we use 
(£ ~ 0.145). The average particle density at the vortex core is 
then given by 



iV(|*o(r C orc)| 2 ) = 



(£ + 2)C( 



— V(-i)" ir 
o ^ v ; 



xe 

xe^(fe) r "- (ai+a2 - 2r - ) 
xe - 2^K+2) a -( r " x ( ai+a2 - 2r ™» 

(ai+a 2 -2r m ) 2 



xe 2! «<«+ 2 > 



(B2) 
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where n = N/A is the mean particle density. 



Appendix C 



Appendix A 

Here we outline how to express the matrix elements in 
terms of analytic functions. We can sum over the first pair of 
variables in Eq.[lO] say mn an d m i2, using the identity rf26tl 



E(- r 



mn^—7rt(m -fn ) — irtmn-\-2mix~\-2niy 



ee d 3 (x\it)ti 3 (2y - x\3it) - i} 1 (x\it)'d 1 (2y - x\3it), 

(Al) 

where d.- t are the Jacobi Theta-functions and t = 1/ How- 
ever, x and y will themselves be dependent on other summa- 
tion variables (m 3 i and m 3 2 in this case), which need to be 
extracted out of the Theta-function arguments before they too 
can be summed over using Eq. lAll They can only be extracted 
if their coefficients are multiples of ir/2 or it/ 2, where it is 
the period of a given Theta-Function. A number of identities 
that exploit the quasi-periodicity of the Theta-functions were 
used in this 'extraction' (see Ref. ll27ll for an exhaustive list). 



We obtain the expectation value of Hi for the trial wave- 
function given in Eq.|28] Let 



U 



N 



exp 



E T K 6 -* 



ke 



(CI) 



and | Atrial) = Un\N), where Un is a depletion operator that 
conserves total particle number. We have made use of the 
following identities which, with the exception of the last one, 
mirror those given in Eq. [26] 

U N b n U N ee 6 q cosh6» q - 6t jSe^" sinh6> q , 



U N b\U N ee b\ cosh6i q - fe-q/r 1 ^ 11 ^ sinh I 



u f N b Q b u N = p^/W^iVn 

+0 (D 2 /N) , 



2N- 1 
2N(N- 1) 



(C2) 



Appendix B 

Here we derive the relation between the vortex core density 
and the mean-squared vortex core displacement. We take the 



where D = X) q ^o bqbq is the depletion operator. Unlike U, 

Un does not commute with &o and b\, and therefore a third 
identity was required. After a somewhat lengthy calculation 
one finds 
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Atrial \Hl Atrial) = Af(0, 0,0,0) 



JV 2 - TV + (1 - 27V) ^2 sinh2 + H sinh2 O + 2 E sinh2 ^q cosh2 ^ 



+4 M (°, q> 0, q) sinh 2 q iV - ^ sinh 2 (9 k - 2 cosh 2 (9 q I + ^ 2 sinh 2 (9 k sinh 2 q A/(k, 0, -k + q, k + q) 

q#0 V k^O / k#0,q#0 



- i \/N 2 -N ^ sinh 26> q 



q#0 



2iV ~ 1 ; | sinh 2 6„ + 1 Y sinh 2 k | - 1 
iV(iV-l)^ q 2^ k J 



- M (q, -q, -q, -2q) + e~^M (0, q, -q, 



^2 sinh6» q cosh q sinh 6» k cosh6> k e i(0k ~ 0c,) M(k, -k + q, -q- k, -2k) + O 

k ^0,q^0 



^Esinh 2 ^ 



r 



When finding the set of parameters 6> k and k that min- other words, we are constrained to optimise in the region of 
imises the expectation value of Hj, we need to check that the parameter space where the incurred error is small, 
error ~ O ((N — Nq) 2 /N) is small and can be neglected. In 
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